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Abstract 

We optimize Hockney and Eastwood's Particle-Particle Particle-Mesh (P3M) algorithm to 
achieve maximal accuracy in the electrostatic energies (instead of forces) in 3D periodic charged 
systems. Fo this end we construct an optimal influence function that minimizes the RMS errors 
in the energies. As a by-product we derive a new real-space cut-off correction term, give a trans- 
parent derivation of the systematic errors in terms of Madelung energies, and provide an accurate 
analytical estimate for the RMS error of the energies. Fhis error estimate is a useful indicator 
of the accuracy of the computed energies, and allows an easy and precise determination of the 
optimal values of the various parameters in the algorithm (Ewald splitting parameter, mesh size 
and charge assignment order). 



I. INTRODUCTION 



Long range interactions are ubiquitously present in our daily life. The calculation of 
these interactions is, however, not an easy task to perform. One needs indeed to resort to 
specialized algorithms to overcome the quadratic scaling with the number of particles, as 
soon as the simulated system includes more than a few hundred particles, see for example the 
review of Arnold and HolnP. In Molecular Dynamics simulations, one is mainly interested 
in the accuracy of the force computation, since they govern the dynamics of the system. In 
contrast, in Monte Carlo (MC) simulations, the concern is to compute accurate energies. 
If the potential is of long range (e.g. a Coulomb potential or dipolar interaction), and 
one has chosen to use periodic boundary conditions, the computation of both observables 
is quite time consuming if one uses the traditional Ewald sum. Since the seminal work 
of Hockney and EastwoocP it has been common to resort to a faster way of calculating 
the reciprocal space sum in the Ewald method with the help of Fast-Fourier- Transforms 
(FFTs). These algorithms are called mesh-based Ewald sums, and various variants existP. 
They all scale as iVlogiV with the number of charged particles N, and the algorithms are 
nowadays routinely used in simulations of bio-systems, charged soft matter, plasmas, and 
many more areas. The most accurate variant is still the original method of Hockney and 
Eastwood, which they called particle-particle-particle-mesh (P3M), and into which various 
other improvements like the analytical differentiation used in other variants of the mesh- 
based Ewald suirP can be built in. In addition, an accurate error estimate for P3M exists, 
so that one can tune the algorithm to a preset accuracy, thus maximizing the computational 
efficiency before doing any simulations^. 

While in the standard P3M algorithrrP, the lattice Green function, called the "influence 
function", is optimized to give the best possible accuracy in the forces, the electrostatic 
energy is usually calculated with the same force-optimized influence function. However, 
there are certainly situations where one needs a high precision of the energies, for instance 
in Monte Carlo simulations, and the natural question arises whether one can optimize the 
influence function to enhance the accuracy of the P3M energies. The main goal of this paper 
is to derive the energy-optimized influence function, and to derive an analytical estimate for 
the error in the P3M energies. This error estimate is a valuable indicator of the accuracy 
of the calculations and allows a straightforward and precise determination of the optimal 
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values of the various parameters in the algorithm (Ewald splitting parameter, mesh size, 
charge assignment order). 

The present derivation of the optimal influence function, and the associated error esti- 
mate, is concise and entirely self-contained. The present paper can thus also serve as a 
pedagogical introduction to the main ideas and mathematics of the P3M algorithm. 

The paper is organized as follows. In Sec. [TTJ we briefly review the ideas of the standard 



Ewald method and provide the most important formulae. In Sec. |III[ we derive direct 
and reciprocal space correction terms which compensate, on average, the effects of cut-off 
errors in the standard Ewald method. We interpret the formulae in terms of the direct and 



reciprocal space components of the Madelung energies of the ions. In Sec. [TV] the calculation 
of the reciprocal energy according to the P3M algorithm (i.e. with a fast Fourier transform 
and an optimized influence function) is presented. The mathematical analysis of the errors 
introduced by the discretization on a grid is performed in Sec. |V| This analysis is used in 



Sec. VI to derive the energy-optimized influence function and the associated RMS error 
estimate. The derivation shows that the P3M energies must be shifted to compensate for 
systematic cut-off and aliasing errors in the Madelung energies of the ions. Finally, our 



analytical results are tested numerically in Sec. VII 



II. THE EWALD SUM 

We consider a system of N particles with charges qi at positions rj in an overall neutral 
and (for simplicity) cubic simulation box of length L and volume V = L 3 . If periodic 
boundary conditions are applied, the total electrostatic energy of the box is given by 

N 

^=2EE wK + nL), (2.1) 

where v(r) = l/\r\ is the Coulomb potential, ry = r« — rj, and n is a vector with integer 
components that indexes the periodic images. The prime indicates that the (divergent) 
summand for i = j has to be omitted when n — 0. 



Because of the slow decay of the Coulomb interaction, the sum in (2.1 ) is only condition- 
ally convergent: its value is not well defined unless one specifies the precise way in which 
the cluster of simulation boxes is supposed to fill M. 3 . Often, one chooses a spherical order 
of summation, which is equivalent to the limit of a large, spherically bounded, regular grid 



of replicas of the simulation box, embedded in vacuum. The simulation box can then be 
pictured as the central LEGO brick in a huge ball made up of such bricks. If this "lego ball" 
is surrounded by a homogeneous medium with dielectic constant e' (e' = 1 if it's vacuum) 
and if the simulation box has a net dipole moment M = ^2 i qiVi, the particles in the ball 
will feel a depolarizing field created by charges that appear on the surface of the uniformily 
polarized ball. It can be shown that the work done against this depolarizing field when 
charging up the system is 

&*> = , 2? ™ 2 ■ (2.2) 
(1 + 2c')L 3 y 1 

in the case of a spherical order of summation^ (for other summation orders, see the articles 

of Smith 8, and Ballenegger and HanserP). The energy is contained, even if not easily 



seen, in the total electrostatic energy (2.1) (at least when e' = 1 since such a vacuum 



boundary condition was assumed in writing <\2.1 )). Obviously, the energy E^ vanishes if 



we employ metallic boundary conditions defined by e' = oo. 

The fact that E^ depends on the order of summation, and hence on the shape of macro- 
scopic sample under consideration, is a consequence of the conditional convergence of the 



sum (2.1). Due to the energy cost E^ d \ the fluctuations of the total dipole moment of the 
simulation box (and hence of the considered macroscopic sample) depend on the dielectric 
constant e' and on the shape of the sample. The energy E^ is crucial to ensure, for ex- 
ample, that the dielectric constant e of the simulated system obtained from the Kirkwood 
formula^, which relates e to the fluctuations of the total dipole moment, is independent of 
the choices made for the sample shape and for the dielectric boundary condition 1 11 1 12 1 . 



Ewald's method to compute the energy (2.1 ) is based on a decomposition of the Coulomb 



potential, v(r) = ip{ r ) + 0( r )> such that ip(r) contains the short-distance behavior of the 
interaction, while 4>(r) contains the long-distance part of the interaction and is regular at 
the origin. The traditional way to perform this splitting is to define 

4>(r) = erf (ar)/r, r = \r\, (2.3) 

and 

ip(r) = v(r) — <f)(r) = erfc(ar)/r (2.4) 

With this choice, ip(r) corresponds to the interaction energy between a unit charge at a 
distance r from another unit charge that is screened by a neutralizing Gaussian charge dis- 



tribution whose width is controlled by the Ewald length a 1 . Following this decomposition 
of the potential, the electrostatic energy can be written in the well-known Ewald formP^ 

E = E {r) + E (k) + E (d) (2.5) 

where the real-space energy E^ contains the contributions from short-range interac- 
tions i[>(r), i.e. 

N 

n&? i,j=l 

and the reciprocal space energy E^ k ' contains contributions from long-range interactions 0(r) 
(apart from the contributions that are responsible for the conditional convergence which are 
included in the term E^ in (2.5)). The fact that the surface term (or "dipole term") E^ 



is independent of the Ewald parameter a shows that this contribution is not specific to 
the Ewald method, but more generally reflects the problems inherent to the conditional 



convergence of the n sum in Eq. (2.1). Contrary to E^ which can be computed easily in 
real space thanks to the rapid decay of the ip interaction, E^ k ' is best computed in Fourier 
space, where it can be expressed as^ 

E (k) = E (ks) _ E {s) 

where 

fceK 
fc^o 

= Q 2 ^- (2.9) 

\/7T 



with 

N 

Q 2 = ^l ( 2 - 10 ) 



1=1 



In (2.8), 4>(k) is the Fourier transform of the reciprocal interaction (2.3), 



0(fc) = J e - lfc ' r 0(r)dr = ^exp(-A; 2 /4a 2 ), (2.11) 
and p(k) is the Fourier transformed charge density 

N 

p(fc) = $>e-^. (2.12) 

i=i 



The sum in (2.8) is over wave vectors in the discrete set K = {2-7rn/L : n e Z 3 }. The term 
k = is excluded in the sum because of the overall charge neutrality. The self-energy term 
E^ compensates for the self-energies (the reciprocal interaction of each particle with itself 



lqf(j){r = 0) = qfa/y/n) that are included in E^. 



The energy (2.1) converges only for systems that are globally neutral. For systems with 
a net charge, the sum can be made convergent by adding a homogeneously distributed 
background charge which restores neutrality. In that case, an additional contribution^ 

N 



1=1 



must be added to (2.5) to account for the interaction energies of the charges with the 
neutralizing background. 

The reciprocal energy E^ ks \ defined by the Ewald formula (2.7), is the starting point of 



mesh-based Ewald sums, which are methods to compute efficiently that energy in many- 



particle systems. Notice that (2.7) can also be written in an alternative form in terms of a 
pair potential and a Madelung self-energy, see Appendix |A} The inverse length a tunes the 
relative weight of the real space E^ and the reciprocal space E^ k > contributions to the energy, 
but the final result is independent of a. In practice, E^ and E^ can be computed using 



cut-offs, because the sum over n in (2.6) and the sum over k in (2.8) converge exponentially 



fast. Typically, one chooses a large enough to employ the minimum image conventional in 



Eq. fl2.6[ ). 

At given real and reciprocal space cut-offs r cut and k cnt , there exists actually an optimal a 
such that the accuracy of the approximated Ewald sum is as high as possible. This optimal 
value can be determined with the help of the estimates for the cut-off errors derived by 
Kolafa and PerramP^, by demanding that the real and reciprocal space contributions to the 
error are equal. Kolafa and Perram's root-mean-square error estimates are 

-« 2 r?„ t 



^^l^hk^ (2 ' 14) 



and 



e -(nk cut /aL) 2 



A£^Q 2 « ■ (2.15) 

^ ^cut 

These error estimates make explicit the exponential dependence of the error on the real and 
reciprocal space cut-offs. 



Formula (2.15) is actually valid only when a correction term (given by Eq. (3.8) below) 
is added, to compensate the systematic error that affects the reciprocal energies when the 



sum over wave vectors in (2.8) is truncated. The origin of this correction term is explained 



in detail in Sec. Ill A similar term must also be introduced in the P3M algorithm when 



one computes the electrostatic energy. Similarly, the direct-space energy (2.6) also contains 
a systematic error when the pair- wise interaction is truncated at the cut-off distance r cut . 
The derivation in the next section will also provide a correction term for this effect. 
To summarize, the final Ewald formula for the total electrostatic energy reads 



E = E ir) 


[eq.(2.6) 




[eq.(2.8) 


_#(.) 


[eq.(2.9) 


+ E W 


[eq.(2.2) 


+ E {n) 


[eq.(2.13) 



(2.16) 



Furthermore, when the sums in E^ and are evaluated numerically using cut-offs, an 



additional correction term E cut , defined in Eq. (3.13) below, must be added to the truncated 
energy, as shown in the next section. 



III. CORRECTION TERM FOR TRUNCATED EWALD SUMS 

If we consider electroneutral systems where the charged particles are located at random, 
we expect the electrostatic energy to vanish on average, because there is an equal probability 
to find a positive or negative charge at any relative distance r. However, when periodic 
boundary conditions (PBC) are applied, the average energy of random systems does not 
vanish, because each charge interacts with its own periodic images (and with the uniform 
neutralizing background provided by the other charges). 

Since this interaction energy Ei mg of an ion with its periodic images and with the neu- 
tralizing background does not depend on the position of the ion in the simulation box, it 
plays the role of a "self-energy". We will refer to E img as the Madelung (self-)energy of an 
ion, to avoid confusion with the self-energy |q ,2 0(O) already defined in the Ewald method as 
the reciprocal interaction of a particle with itself. 

We denote by angular brackets the average over the positions of the N charged particles: 
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A. Madelung energy 

The Madelung energy of an ion takes the form E{ mg = \q 2 Ci where C, is a purely numerical 
factor in units oil/L that depends only on the size and shape of the simulation box. 

Let us calculate the average electrostatic energy of random charged systems in PBC, 
to find the value of ( and derive a correction term for cut-off errors in truncated Ewald 



sums (some results derived here will be used in Sec. VIA). On the one hand, the average 
Coulomb energy of the random systems is by definition Q 2 (/2, while on the other hand, it 
can be calculated as the sum of a direct space contribution (E( r >} and a reciprocal space 



contribution (E^ k '\. The average reciprocal energy is, using (2.9) and (2.8), 



2L 3 — - — - ' WTT 

i,j keK v 

fc^O 

Since (exp(— ik ■ rj)) = Sk.o, all terms with j ^ i vanish (this is due to the fact that the 
Ewald pair potential averages to zero, see Appendix [AJ . By contrast, "self" terms (i = j) 
remain and lead to 

<* w > = f (^E*M--%) = t< w . (3 ' 3) 

fc^O 

where the second equality defines The average real-space energy of a single ion of 

charge qi in periodic random systems is 

E\ r) ) = f ( £ V>(nL) - ^ / ^(r)d 3 r) (3.4) 

where the first term is the sum of the direct interactions of the ion with all its periodic 
images, while the second term corresponds to its interaction with the uniform background 
charge density —qi/L 3 provided by the other particles in the system. Since 

/ ip(r) d 3 r = Air [ r 2 tp(ar) dr = (3.5) 
Jr3 Jo a 
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we can write the average total real-space energy as 

<£ M > = ?(£*("£>--^) = ?C<", (3.6) 



2 \ ^ rv ' a 2 L 3 J 2 
which defines C ■ The second term in is, not surprisingly, identical to the energy 



defined in (2.13). Notice that the above result for (^E^ r '\ may also be obtained by splitting 
(2.6) into self (i = j) and interaction terms, and using for the latter Ylj^iQj = — <?« which 
follows from the electro-neutrality condition. The expression of the factor £ = Q r ' + is 
therefore 



Eq. (3.7) can be computed for a number of different box geometries^. For a cubic 



simulation box of size L, it yields^^ 

C ~ -2.837297479480619610825442578061/L. 
The above calculation shows that, when a charged system is simulated using PBC, the 



electrostatic energy (2.1) contains an additional contribution Q 2 (/2. The existence of this 
Madelung self-energy can be made more apparent in the Ewald formula for E, as shown in 
Appendix [A} 

B. Madelung cut-off error correction terms 



The Ewald sums (2.6) and (2.8) are necessarily truncated when evaluated in a simula- 
tion. These truncations introduce systematic cut-off errors in the total energy, because the 
Madelung self-energies of the ions are then not fully accounted for. This systematic error is 
typically of the same order of magnitude, or even larger, than the fluctuating error, due to 
the use of cut-offs, in the Ewald pair interaction energy^SEfi]. Note, that no similar systematic 
error affects the electrostatic forces, because the Madelung energy does not depend on the 
position of the ion. 

Fortunately, it is easy to suppress the systematic bias in the computed energies. We 
simply have to add the cut-off correction 

fceK 

k^>hcut 



to the computed fc-space energies, which Kolafa and Perram termed the diagonal correc- 
tion^. The value of E^ t does not depend on the configuration and may thus be computed 



in advance using a sufficiently large second cut-off k' cut > k > k cut . Using definition (3.3) 
of C {k) , 

we can rewrite (|3.8l) as 



E (k) 

^cut 



o_ 

2 



m - c {k) 

S Scut 



(3.9) 



where 



Ak) 

*5CUt 



2a 



h - 5§- (3 - 10) 

fe^O, k<k cut 

Similarly, if the real-space energies are computed using a cut-off r cut < L/2 (minimum 
image convention), we see from Eqs. (3.4), (3.5), and (3.6), that the r-space cut-off correction 



^ul = ?(c (r) -d u b (3.ii) 



2 

where 



Wo 



Tout 



Cut = / r ^(r)dr 



must be applied to the direct space energies. It is natural that the correction terms E cn [ 
and E^ u \ are made up of the exact Madelung energies, minus the average Madelung energies 
of the ions as obtained from a calculation with direct and reciprocal space cut-offs r cut and 

^cut- 



Adding (3.9) to (3.11) and using (3.7), the two cut-off corrections can be combined into 



a single expression 

-C'cut — -C'cut + ^cut — — Scut ~ Scut ) • [6.16) 



2 

All of these terms can easily be precomputed numerically before the start of a simulation. 
Correcting the systematic cut-off errors in the energies with the term E cut does improve 
significantly the accuracy of the results, especially when working with small cut-offs. In 

(r) 

numerical tests, however, the direct space cut-off correction E^. ut has been found to be 
mostly negligible compared to the reciprocal space correction E^ t for all practical purposes. 
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IV. MESH-BASED EWALD SUM 



The idea of particle-mesh algorithms is to speed up the calculation of the reciprocal 
energy E^ ks ' with the help of a Fast-Fourier- Transform (FFT). To use a FFT, the charge 
density must be assigned to points on a regular grid. There are several ways of discretizing 
the charge density on a grid, and to get the electrostatic energy from the Fourier transformed 
grid. We will use the P3M method of Hockney and Eastwood (but with the standard Ewald 
reciprocal interaction ( 2.3[ )), because this method surpasses in efficiency the other variants 



of mesh based Ewald sums (PME, SPME)E1. 

For simplicity, we assume the number of grid points M to be identical in all three direc- 
tions. Let h = L/M be the spacing between two adjacent grid points. We denote by M the 
set of all grid points: M = {mh : m G Z 3 and < rn Xj y iZ < M}. 

The mesh based calculation of the reciprocal energy is made in the following steps: 

A. Assign charges to grid points 

The charge density PM( r ) &t a grid point r is computed via the equation 

p M (r) = Ju(r- r')p{r')dr', r e M, (4.1) 

where U(r) = h~ 3 W(r) with W the charge assignment function (the factor h~ 3 ensures 
merely that PM( r ) has the dimensions of a density). A charge assignement function is clas- 
sified according to its order P, i.e. between how many grid points per coordinate direction 
each charge is distributed. Typically, one chooses a cardinal B-spline for W, which is a 
piece-wise polynomial function of weight one. The order P gives the number of sections in 
the function. In P3M, we only need the Fourier transform of the cardinal B-splines, which 
are 

w(P)(i,\ -h*( sin (W 2 ) sin(fc y h/2) sin(fc 2 V2) ^ P 

W [ >~ V k x h/2 k y h/2 k z h/2 J ■ [ ^ Z) 

Notice that PM( r ) = J2i Qi W(r — Tj), apart at the boundaries where the periodicity 
has to be properly taken into account. 
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B. Fourier transform the charge grid 



Compute the finite Fourier transform of the mesh-based charge density (using the FFT 
algorithm) 



pm 



(fc) = h 3 J2 PM(r)e-^ r = FFT{p M }, fc £ M. (4.3) 



Here fc is a wave vector in the reciprocal mesh M = {2irn/L : n £ Z 3 , |Tia. tJ/)JS | < Af/2}. 

We stress that Pm(^) differs from p{k) for fc £ M, because sampling of the charge density 
on a grid introduces errors (see Sec. M). 



C. Solve Poisson equation (in Fourier space) 

The mesh-based electrostatic potential $m is given by the Poisson equation, which re- 
duces to a simple multiplication in k-space: 

$ M (k) = PM(kMk), fc£M, (4.4) 



with 4>(k) the Fourier transformed reciprocal interaction (2.11). However, instead of using 



4>(k) in the above equation, it is better to introduce an "influence" function G(k). We 



replace therefore Eq. (4.4) by 



^ M (k)=p M (k)G(k), fc£M. (4.5) 

where G(k) is determined by the condition that it leads to the smallest possible errors in 
the computed energies (on average for uncorrelated random charge distributions). G(k) will 



be determined later (see Eq. (6.21)); it can be computed once and for all at the beginning 
of a simulation since it depends only on the mesh size and the charge assignment function. 
G{k) plays basically the same role as the reciprocal interaction 0(fc), except that it is tuned 
to minimize a well defined error functional in pj^(k). We stress that G(k) is defined only 
for fc £ M (we dropped the subscript M on the influence function to alleviate the notation). 
The idea of optimizing G(k), which is a key-point of the P3M algorithm, ensures that the 
mesh based calculation of the reciprocal energy gives the best possible results^ 
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D. Get total reciprocal electrostatic energy 



Expression (2.8) is approximated on the mesh by 



(4.6) 



fceM 
fc^o 



The total reciprocal energy follows from subtracting the self-energies from the above quan- 
titv . =;(*) _ E (ks) _ E ( S ) 

my. -C/ P3M — -C/ P3M ^ ■ 



E. Electrostatic energy of individual charges (optional) 

If the reciprocal energy of each individual particle is needed (and not only their sum as 
in step 4), the potential mesh must be transformed back to real space via an inverse FFT, 
i.e. 

®M(r m ) = ^ E ^M(fc)e' fc - rm = FFT- 1 ^}. (4.7) 
fceM 

The mesh-based potential is then mapped back to the particle positions using the same 
charge assignment function: 

$(r):= W{r-r m )$ M (r m ). (4.8) 

r m SM p 

In this equation, M p = {mh : m G Z 3 } is the mesh extended by periodicity to all space, 



and $M( r ) is assumed to be periodic (with period L). The interpretation of Eq. (4.8) is 
the following: due to the discretization each particle is replaced by several "sub-particles" 
which are located at the surrounding mesh points and carry the fraction W(r — r m ) of the 
charge of the original particle. The potential at the position of the original particle is given 
by the sum of the charge fraction times the potential at each mesh points. The reciprocal 
electrostatic energy of the i th particle is then gj$(r.;)/2, and the total reciprocal energy 
(including self-energies) is the sum 



-^psm 2 



p(r)$(r)dr = i^«*(r 4 ). (4.9) 



This formula gives the same result for the total energy as Eq. (4.6). A mathematical proof 
of the equivalence is given in Appendix [Cj 
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V. ANALYSIS OF DISCRETIZATION ERRORS 



If the fast Fourier transform has the benefit of speed, it has the drawback of introducing 
errors in the /c-space spectrum of the charge density: p~yi(k) differs from the true Fourier 



transform (2.12)(times a trivial factor U(k)) because of the discretization on a finite grid. 

The difference is two-fold. Firstly, p(k) is defined for any vector in the full k-space K, 
whereas pjy[(fc) i s defined only for k e M, i.e. in the first Brillouin zone. This is a first 
natural consequence of discretization: if the grid spacing is h, it necessarily introduces a 
cut-off |fca!,j,,«| < 7T Jh in /c-space. Secondly, the act of sampling the charge density at grid 



points, which is mathematically embodied in Eq. (4.3) by the presence of a discrete Fourier 



transform instead of a continuous FT, introduces aliasing errors. While a continuous FT 



would simply transform the convolution Eq. (4.1) into 



FT{p M }(fe) = U(k)p(k), kEK, (5.1) 
the finite Fourier transform results in (see proof in Appendix |bT) 



p M (k) = FFT{p M }(fc) = J2 U ( k + rnk g )p{k + mk g ), keM, (5.2) 

where k g = 2n/h. The sum over m shows that spurious contributions from high frequencies 
of the full spectrum U (k)p(k) are introduced into the first Brillouin zone M. These unwanted 
copies of the other Brillouin zones into the first one are known as aliasing errors^. 

To avoid aliasing errors, the spectrum needs to be entirely contained within the first 
Brillouin zone. Since p(k) may contain arbitrary high frequencies, this can only be achieved 
by choosing U (k) to be a low-pass filter satisfying U (k) = for k e K \ M. But the 
charge assignment function would then have a compact support in /c-space, and hence an 
infinite support in r-space. This is not acceptable, as it would require the grid to have an 
infinite extension. The need to keep the charge assignment function local in r-space means 
that U(k) cannot be a perfect low pass filter. Aliasing errors are therefore unavoidable, 
and the impact of these errors must be minimized, by choosing a good compromise for the 
charge assignment function and optimizing the influence function. The influence function 
can indeed compensate partially for the aliasing errors, because the spectrum of U(k) is 
known exactly at all frequencies. 

The error in reciprocal energy, for a given configuration p(r) of the charges, is defined by 
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the difference 



AE (k) = E 



(fc) 

P3M 



(5.3) 



where is the exact reciprocal energy (see (2.8) and (2.9)). The above analysis of dis- 



cretization errors results in the explicit formula for this error 

AE(k) = ^ E \pM(k)\ 2 d(k) -^- 3 J2 \p( k )\ 2 m, 



(5.4) 



fceM 
fc^o 



fceK 
fc^o 



where PM.(k) is given by (5.2). The error AE^ is due to the finite resolution h offered by the 
mesh. The finiteness of h introduces the cut-off ir/h in fc-space (k e M) and causes aliasing 
errors (pm(^) 7^ P~{k)U(k)) that cannot be entirely eliminated by the charge assignment 
function. 



VI. OPTIMIZATION OF THE P3M ALGORITHM 



We derive in this section the influence function G(k) that minimizes the error (5.4) on 



average for uncorrelated systems, and give a formula for the associated RMS errors. The 



average over random systems is denoted by angular brackets, as in Sec. [ITT 

Notice that the assumption of the absence of correlations is never satisfied in practice 
(even for uniform systems because negative charges tend to cluster around positive charges 
and vice- versa). The error estimate proves however to predict quite accurately the error in 
real systems with correlations, notably in liquids where the pair distribution function g(r) 
decays rapidly to one. 

A. Shift in the energies to avoid systematic errors 



The P3M energies (4.6) contain in general systematic errors, i.e. (AE^) ^ 0, because 
the Madelung energies of the ions obtained in the mesh calculation contain cut-off and 
aliasing errors. The average error 

K = (AE^) = (E^ M ) - (E^) (6.1) 

is a constant that must be subtracted from the P3M energies, to ensure that the energies 
are right on average. The corrected P3M energies are thus obtained by applying a constant 
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shift to the original P3M energies: 



E (k) 



E {k) - K 



(6.2) 



where the constant K depends on the various P3M parameters like mesh size, charge as- 
signment order (CAO) and Ewald splitting parameter. 

Let us determine analytically the constant (6.1). Writing it as K = (E p k ^ M \ — (E^}, 



we can use the result Q for (#<*)): it is nothing but Q 2 ( (k) /2, i.e. the A;-space Madelung 



energies of the ions. The other term (Ep^ M \ can be calculated in the same way as (3.2) 



/ 



Using (4.6), (5.2) and (2.9), we find 



kelL 



m6Z 3 



2a 



2 ^P3M' 



(6.3) 



(k) 

which defines Cp3m- The result (6.3) can be interpreted as the average fc-space Madelung 
energies of the ions as obtained from the mesh calculation, i.e. including cut-off and aliasing 



errors. The explicit expression of the correction constant (|6.1|) is thus 
Q 2 



K 



Ak) 

SP3M 



c 



(fc) 



Q 2 



2L 3 



(6.4) 



fceM 
fc^o 



kew, 

fc^O 



In the last sum in (6.4), the terms with \k Xj y :Z \ > n/h are equivalent to the fc-space cut- 



off correction defined in (3.8). These terms compensate for the fact that the Madelung 



energies of the ions are underestimated in the mesh calculation because of the cut-off n/h 



introduced by the finite size of the mesh. The remaining terms in (6.4) compensate, on 



average, the aliasing errors that affect the Madelung energies of the ions obtained from the 
mesh calculation. 

(r) 

Notice that the two correction terms E^ t and — K can be combined together in the simple 
expression 

Q 2 



rpcut _ rp( r ) _ U — II_ ( f _ A k ) _ A r )\ 
^VZM — -^cut 2 V ^P3M Scut J 



(6.5) 



— (r) 

where ( is defined by (3.7) and Q u [ is given in (3.12). We stress that E- 



7P3* M has the same 



structure as the correction term (3.13) for truncated Ewald sums. The difference lies in the 



replacement of Qut by the quantity Cp3m defined in (6.3), which accounts for both the cut-off 



and aliasing errors that affect the reciprocal energies computed on the mesh. 
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In summary, the final formula for computing the total electrostatic energy with the P3M 
algorithm is 

E^E {r) [eq.ra 



I J?{ks) 

-^PSM 

_£<«) 



[eq.(4.6)] 



[eq.(2.9|] 



[eq.(2.2) 



(6.6) 



[eq.(2.13) 



+ E" ut 



v P3M [eq.gjp; 

The correction term E^ M is necessary to compensate on average for systematic errors in the 
mesh calculation. It can be computed once for all before the start of a simulation, since it 
depends only on the size of the simulation box, the size of a mesh cell, the charge assignment 
function and the influence function. 



B. RMS error estimate for energy 



The result (5.4) is an exact measure of the error in the P3M energies for a given configu- 



ration p(r) of the particles. Let us average this expression over all possible positions of the 
particles to get a useful overall measure of the accuracy of the algorithm. The RMS error 
of the corrected P3M energies is, by definition, 



v rms; 



- E^f) = <(A£7« - Kf) 



(6.7) 



where we used (5.3) and (6.2). We can isolate in AE^ "interaction" terms (i ^ j) from self 
terms (i = j): 

2 X 



(AE 



(*) ^2 
RMS/ 



( AE^l + AE^ - K 



(6.1 



We recall from Sec III that the interaction terms vanish on average for random systems: 



0. The correlation 



(6.9) 



vanishes as well for the same reason (this is due to the fact that the average Ewald interaction 



energy between a fixed particle i and a particle j ^ % is zero, see Appendix A). Eq. (6.8) 
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reduces therefore to 

(A4^s) 2 = ((A4?) 2 ) + ((A£<S " Kf) , (6.10) 

where the first term accounts for fluctuating errors in the interactions energies, and the 
second term accounts for fluctuating errors in the corrected Madelung self-energies of the 
ions. Since the latter term may be written as 

((A^S - Kf) = ((AE^f) - K\ (6.11) 

we remark that the shift —K derived in the previous section, in addition to removing the 
systematic bias in the fc-space energies, also reduces the fluctuating errors of the /c-space 
self-energies by an amount —K 2 . 



In the substraction AE^ f — K, it can be seen, from (5.4) in which only i — j terms are 



kept and (6.4), that all terms containing <p(k) cancel out, so we have 



fceM * ™i m 2 

fc^o 



-;>> 2 (fc™)l) ) (6.i2) 



where we used the symmetry U(—k) = U(k) and introduced the shorthand notation k m = 
k + k g m. When the square is expanded, the summation over particles Y2i becomes a double 
summation Yin'- AH terms with i' ^ i vanish, because (exp(ik g (m 1 — m 2 ) ■ 7*j)) = S rni>m2 , 
leaving identical sums over m which cancel each other. The remaining terms i' = i evaluate 
to 



mi 7TO2 T713 



with 



i fceMfc'eM 
fc^o feVo 

m m! 

= i^(E^self (6-13) 

i 

^ 2 eif = ^E E^^'HE E E^^)^^)^^)^-!--^)}- 

fceMfc'eM ™i m 2 ^mi ra 3 

fc^O fcVo 
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(6.14) 

The fluctuating errors of the Madelung self-energies scale therefore like Yli it with the 
valencies of the ions. The prefactor is somewhat complicated since it involves a double 
summation over wave vectors and a triple summation over alias indices mi,m.2,m3, but 
H^ clf can be evaluated reasonably fast. The numerical calculation of H^ clf can be accelerated 
by taking profit of the symmetries (the sum over k can be restricted to only half an octant of 
the reciprocal mesh), and by skipping inner loops in the triple summation over alias indices 
if the product of the charge fractions is almost zero. 

We calculate now the fluctuations of the errors in the interaction energies, i.e. the first 



term of Eq. (6.10). That term reads, using (6.8), (5.4), (5.2) and (2.12) and keeping only 



interaction terms: 



ikm -Tj 



fceM V 

fc^O i¥=3 



X 



G{k) e- ikm " ri U(k m )U(k m ,) - e-^'ftk, 
Tri'ez 3 



. (6.15) 



The calculation of this average is straightforward, though somewhat tedious. We find that 
it reduces to 

Q 4 



(k)\2 
int ) 



4L 3 



H 



int' 



(6.16) 



where 



fceM 

fc^O 



G\k) ( U\k m )) - 2G(k) VfrmMkm) + E <P\k r 



. (6.17) 



The factor 2 in Hf nt originates from the fact that each pair of particles appears twice in 



the sum over % and j ^ i in (6.15). Expression (6.17) is the analog for the energy of the 



parameter Q introduced by Hockney and Eastwood to measure the accuracy of the P3M 
forces^. Notice that (6.17) is given in real space by 
2 



Hi 



hit 



V, 



cel1 JVceU 



dri / dr 



>P3M 



L 3 



V;ri) - (j)(r - n)Y 



(6.18) 



where 0p3m(^;^i) is the reciprocal potential at r created by a unit charge located at ri, as 
obtained from the P3M algorithm. (This potential is given in Fourier space by combining 
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(C4) with (5.2) in which we set p(r) = S(r — Ti).) H? nt is hence twice the squared deviation 
between the potential 0P3M obtained from the mesh calculation and the exact reciprocal 
potential 0, summed over all relative positions r within the simulation box, and averaged 
over all possible positions of charge r± in a mesh cell (V ce n = h 3 ). 



Inserting the above results in (6.10), our final expression for the RMS error of the (cor- 
rected) P3M energies is 



AE, 



(k) 

RMS 



2 

self 



2L 3 / 2 



(6.19) 



where H 2 nt and H 2 el{ are defined in Eqs. (6.17) and (6.14). This error depends on the 
influence function G{k). The optimal influence function (the one that minimizes the error) 
will be determined in the next section. The above error estimate, together with the optimal 



influence function (6.21) and the constant shift (6.4) which must be applied to the P3M 
energies, constitute the main results of this paper. 



The RMS error (6.19) displays two different scalings with the valencies of the ions: 
(YliQi) 2 f° r errors coming from pair interactions (such a scaling also governs errors in P3M 
forces^) and ^ qf for errors in Madelung self-energies. Because of these different scalings, 
the errors from pair interactions are expected to dominate in systems with many charged 
particles (Q 4 3> Ylii^t)- Notice that H 2 eli is, roughly speaking, proportional to (^2, k G{k)) 2 , 
while itself sca l es hke 'Yl lk G 2 {k). The errors in the Madelung self-energies increase there- 
fore more rapidly than the errors in the pair interaction energies when the Ewald splitting 
parameter a (and hence G(k)) is increased, or when the size of the mesh is increased. The 
importance of the two source of errors (fluctuations in pair interaction energies versus fluc- 



tuations in Madelung self-energies) will be compared in Sec. VII for a test system with 
Q 2 = 100. 



C. Optimal influence function 

We can now determine the optimal influence function G(k), by imposing the condition 



that it minimizes the RMS error (6.19). Since the errors coming from pair P3M interactions 
are expected to dominate the self- interaction errors (except in systems with few particles), 
we optimize the influence function only with respect to the pair interactions. Setting 

Stt2 

^J^- = 0, (6.20) 
6G(k) 
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gives immediately 
G(k) 



£ U\K 



(6.21) 



where we recall that the Fourier-transformed reciprocal interaction (f)(k) is given by (2.11). 
An optimization of the influence function with respect to the full RMS error could be 
performed, but would require solving a linear system of M 3 equations to compute G(k). 



The numerical results shown in Sec. VII will confirm that such a full optimization is not 
necessary in typical systems. 

Since <fi{k) decays exponentially fast, the optimal influence function is given in good 
approximation by 

U 2 (k) 



G(k) ~ <p(k) 



Em,eZ 3 U 2 (k, 



(6.22) 



G(k) differs thus from <p{k) by a factor which is always less than one. This damping of the 
interaction compensates as well as possible for the aliasing errors introduced by the use of 
a fast Fourier transform. If U(k) were a perfect low-pass filter (U(k m ) = if m ^ 0), no 
aliasing error would occur and the influence function would reduce to G(k) = <p(k) / U 2 (k) . 



This is indeed the result expected from (5.4) when aliasing errors are absent. The true 



optimal influence function (6.21) differs from this simple expression by contributions from 



the high-frequency spectrum of U(k) and reciprocal interaction (2.11). 



Hockney and Eastwood obtained the following optimal influence function by minimizing 
the errors in the forces instead of the energy^ 

Em(fc-fcm) U 2 {k m )4>{k m ) 



G {iovces \k) 



k 2 (E m U 2 (k m ) 



(6.23) 



Obviously, this function is also given in very good approximation by (6.22). This explains 



why influence functions (6.21 ), (6.22) and (6.23) all give very similar results when computing 
energies and forces. 



Inserting (6.21) into (6.17), we find that the minimal value of H 2 nt is 

2~ 



int 



-y 

L 3 ^ 



fceM 
fc^o 



£ <P\k 7 



E^ 2 (fcm)0(fc, 

E m u 2 (k m ) 



(6.24) 
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This is the expression of Hf Qt to be used in the RMS error estimate (6.19) when the P3M 



algorithm is optimized to yield the smallest possible errors in the pair interaction energies. 

We recall that the errors in the P3M energies originate from aliasing effects (due to the 
sampling on a grid) and truncation errors (due to the fact that the reciprocal mesh contains 
only a finite number of wave vectors). The truncation error can only be reduced by choosing 
a larger mesh or by using a reciprocal interaction with a faster decay in fc-space, whereas 
the aliasing errors may be reduced by increasing the order of the charge assignment function 
(up to the maximum order allowed by the size of the mesh) . The intrinsic truncation error 



of a given mesh and reciprocal interaction can be obtained by assuming U(k) in (6.24) to 
be a perfect low-pass filter: 



H 2 

int,cut— off 



r 3 L^i 



L 3 



fceM \-mez 3 

fc^O 



(6.25) 



By inserting this formula in (6.19), we get an estimate of the intrinsic RMS cut-off error 



in /c-space, caused by the finite number of wave vectors in the reciprocal mesh. The RMS 



error associated with (6.25) depends only on the size of the mesh and on the choice of the 



reciprocal interaction, i.e. Ewald parameter a if the standard form (2.3) is used. 



VII. NUMERICAL CHECK OF ACCURACY 



In this section, we test the analytical results (optimal influence function, energy 
shift -Ep3M, RMS error estimate) derived in the previous section. We do this by comparing 
the P3M energies with the exact energies calculated in a specific random system. In the 
following, all dimensions are given in terms of the arbitrary length unit £ and charge unit 
C. In particular, energies and energy errors are given in units of C 2 /C. We choose the 
same test system as the one defined in Appendix D of Deserno and HolnP: 100 particles 
randomly distributed within a cubic box of length L = 10C, half of them carry a positive, 
the other half a negative unit charge. The statistical average (• • •) is calculated by aver- 
aging over at least 100 different configurations of this test system (these configurations are 
determined by using the same random number generator as in Deserno and HolnP). Well 
converged Ewald sums (in metallic boundary conditions) were used to compute the exact 
energies of the test systems. The first three systems have energies —15.43059, —15.26641 
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and —15.59147 respectively, values that are all quite close to the Madelung energies of the 
ions Q 2 (/2 - -14.187. 

The P3M energies of the test systems were computed with various mesh sizes (M = 
4,8,16,32), a real-space cut-off r cut = 4.95, and different orders of the charge assignment 
function (from 1 to 7). Our calculations show that using the energy-optimized influence 



function (6.21), instead of the force-optimized influence function (6.23), leaves the energies 
almost unchanged. (A slight improvement in accuracy appears only when the aliasing error 
are at their maximum, namely for a charge assignment order of 1 and large values of a.) 
This behavior could have been expected, since both influence functions are almost equivalent 
to the simple formula ( |6.22[ ). 

We compare in Fig. [I] the measured systematic error (-Ep3M — E exact ) of the uncorrected 
P3M energies [Eq. ( 6.6[ ) without term E^f M ], for CAO's ranging from 1 to 7, to the expected 
bias —EpsM- The agreement is perfect for all CAO's and for all values of Ewald's splitting 



parameter a. The energy shift Ep£ M in Eq. (6.6) removes therefore entirely the systematic 
error, as it should. 

Fig. 1 illustrates that the systematic errors in the uncorrected energies, which are due 
to cut-off and aliasing errors in the Madelung self-energies of the ions, have two different 
contributions of opposite sign. At small values of a, the r-space cut-off error dominates and 
leads to an overestimation of the energy because the negative interaction energy of an ion 
with the neutralizing background charge provided by the other particles is not fully taken 



into account. The cut-off correction (3.11) derived in Sec. Ill does compensate very well for 
this effect. At large values of a, fc-space cut-off and aliasing errors dominate, and lead to an 



underestimation of the Madelung self-energies (expression (6.4) is indeed always negative). 

Since the systematic error {Ep^ M — -EgxLt) i n the reciprocal energies arise solely from self 
terms (the Ewald interaction between a pair of particles is zero on average), this error can 
alternatively, and more efficiently, be measured by computing the P3M energy of a system 
made up of a single ion in the box, averaging that energy over different positions of the 



particle relative to the mesh. To restore electro-neutrality, the interaction energy (2.13) 
with the (implicit) neutralizing background must of course be taken into account before 
comparing the result with the exact Madelung self-energy q 2 (/2 of the ion. This method 
allows one to measure very rapidly the reciprocal contribution to the average error in the 
Madelung self-energies of the ions. We stress that the numerical results shown in Fig. 1 can 
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a 

FIG. 1: Comparison between the measured systematic error of the uncorrected P3M energies 
(crosses) and the theoretical prediction — (solid lines) as a function of Ewald parameter a. 

The average is performed over 1000 test systems consisting in 100 charges located at random in a 
box of size L = 10. The mesh has size M = 8 and the real-space cut-off is r cu t = 4.95. 



easily be transposed to any cubic system with an arbitrary number of ions since the energy 
shift Epf u scales merely as Q 2 /L. 

Having validated the energy shift (|6.5l), we test now the accuracy of the RMS error 



estimate (6.19). We show in Fig. |2j the theoretical predictions for the RMS error of the 
corrected P3M energies for different mesh sizes (thick solid lines), at fixed CAO 2. The 
dominant error at small values of a comes from the truncation in the real-space calculation, 
while fc-space cut-off and aliasing errors dominate at large values of a. The plot shows 

(k) 

also separately the contribution Ai£ self , which accounts for fluctuating errors in the fc-space 
Madelung self-energies, and the contribution AE^ which accounts for fluctuating errors in 
the P3M pair interaction energies. Near the optimal value of a, the error AE^ dominates 
slightly Ai? S gj f by half an order of magnitude. This validates the use of the optimal influence 



function (|6.21|), which was designed to minimize errors in the pair P3M interaction energies 
^seif overcomes AE^ 



(k) (k) 

only. Notice that AE^{ { overcomes Ai? int at large values of a, in agreement with the scaling 



with a discussed in Sec. VI The errors in Madelung self-energies must therefore by included 
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to predict correctly the full RMS error curve in our test system with 100 charged particles, 
but they are expected to become negligible when the number of ions is increased above a 
few hundred. 




FIG. 2: Theoretical predictions for the RMS error of the corrected P3M energies for CAO 2 
and three different mesh sizes [thick solid lines], for the same system and real-space cut-off as in 
Fig. [TJ The two contributions which make up the total fe-space error are also shown independently: 



RMS error in the pair P3M interaction energies (Eq. (6.16), thin solid line) and RMS error in the 



Madelung self-energies (Eq. (6.13), dashed line). 



The predicted RMS errors agree very well with the measured RMS errors, as shown in 
Fig. 3. The small deviations at low values of a are due to a loss of accuracy of Kolafa 



and Perram's r-space error estimate (2.14), and to the fact that this error estimate does 
not take into account the improvement in accuracy brought by the new cut-off correction 



term (3.11). In the regime where the dominant error comes from the /c-space calculation, 
the agreement with our RMS error estimate is excellent, especially at high values of the 
charge assignment order. The errors in the fc-space calculation are caused by truncation 
and aliasing effects. The aliasing errors can be reduced by increasing the charge assignment 



order, but the accuracy cannot go below the minimum fc-space cut-off error (6.25) (dashed 
curve in Fig. 3), which is intrinsic to the mesh size and choice of reciprocal interaction. 
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a 



FIG. 3: Comparison between the measured (crosses) and predicted (solid lines) RMS errors of the 
(corrected) P3M energies, for the same test system, mesh size and real-space cut-off as in Fig. 1. 
The minimal error due to direct and reciprocal space cut-offs is shown as a dashed line. 

The pronounced minimum in the RMS error curves stresses the importance of using the 
optimal value of a when performing simulations with the P3M algorithm (or with the other 
variants of mesh based Ewald sums). Our accurate RMS error estimate for the P3M energies 
can be used to quickly find the optimal set of parameters (mesh size, charge assignment 
order, Ewald splitting parameter) that lead to the desired accuracy with a minimum of 
computational effortP. Whatever the chosen parameters, it can serve also as a valuable 
indicator of the accuracy of the P3M energies. 

VIII. CONCLUSIONS 

In this article, we discussed in detail which ingredients are necessary to utilize the P3M 
algorithm to compute accurate Coulomb energies of point charge distributions. The usage 
of a nearly linear scaling method (~ iVIogiV) like P3M is almost compulsory for systems 
containing more than a few thousand charges. 

In particular, we derived the cut-off corrections for the standard Ewald sum transparently 
and interpreted the systematic errors in terms of Madelung energies. This route lead us to 
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an additional real-space cut-off correction term that has so far not been discussed in the 
literature. Building on these results, we have deduced the fc-space cut-off correction term in 
the case of the P3M algorithm, where additional aliasing errors play a role. Furthermore we 
derived the exact form of the influence function that minimizes the RMS errors in the ener- 
gies, and showed that this function is not much different from the force-optimized influence 
function, which a posteriori justifies why in most P3M implementations the usage of the 
force-optimized influence function does not lead to inaccurate results. Based on the energy 
optimized influence function we derive an accurate RMS error estimate for the energy, and 
performed numerical tests on sample configurations that demonstrate the validity of our 
error estimates and the necessity to include our correction terms. We also demonstrated 
that the electrostatic energy of an individual particle in the system can be obtained in the 
P3M method, but at the expense of an additional inverse fast Fourier transform. 

With the help of the newly derived error estimates we can easily tune the desired accu- 
racy of the P3M algorithm and find suitable parameter combinations before running any 
simulation. 

The P3M algorithm can be generalized along our discussed lines to compute other long 
range interactions. Of particular interest are dipolar energies, forces and torques, and the 
associated error estimates for these quantities. This will be the content of a forthcoming 
publication. Our P3M generalization for the energies will be included in a future version of 
the molecular simulation package Espresso^, that is freely available under the GNU general 



public license. The website http : / /www . espresso . mpg . de provides up-to-date information 
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APPENDIX A: EWALD PAIR POTENTIAL AND MADELUNG SELF-ENERGY 

The Ewald formula for the electrostatic energy E of a periodic charged system can be 
written in a form that underlines the fact that E includes the Madelung self-energies Q 2 C s f2 
of the ions [Q 2 = ^i<7i an d C is defined in (3.7)]. We recall from Sec. [n]that the Ewald 
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formula for E reads, if the system is globally neutral and if we employ metallic boundary 
conditions, 

E = \HY, ft^(r« + nL) + ^ Y, E e~ iHn - ri) m - Q 2 ^=. (Al) 

ij neZ 3 id keK v 

fc^O 

The "self-energy terms" in E, i.e. term E^ and terms i = j, are 
We can write therefore 



E = \ E *<& ( E^ + nL ) + is E e- ifc - r ^(fc)) + y (c + 



7T 



a 2 L 3 , 

1 Q 2 

= 2 E ® & ^ Ewald ( r y ) + "y ^ ( A3 ^ 

where we defined the Ewald pair interactiorP 

^Ewaid(r) = $>(r + nL) + ±^ e -«-0(fc) - (A4) 



Notice that in writing (A3), we used ^i^Zj^Qigj ^ / (c^ 2 L 3 ) = —Q 2 ir/(a 2 L 3 ) which follows 
from electro-neutrality. Thanks to the inclusion of this constant in the definition of VEwaid(f), 
the Ewald pair potential does not depend on the parameter a [d/daVE wa id( r ) — 0] and its 
average over the simulation box is zercP^: 



(^Ewaid(r-)) = ^ j d 3 r ^ Ewald (r) = 0. (A5) 



The latter property is simply a consequence of (exp(ifc • r^)) = 5k,o and Eq. (3.5). 



In conclusion, expression (A3) shows explicitly that the electrostatic energy of a periodic 
charged system includes the Madelung self-energies Q 2 C/2 of the ionsPIS] The fact that the 
Ewald interaction between a pair of particles averages to zero when one particle explores 
the whole simulation box is also noteworthy aspect of Ewald potential. 



APPENDIX B: PROOF OF EQ. (5.2) 



Eq. (5.2 ) is a consequence of the Sampling theorem [refs] and is straightforward to demon- 



strate. The sum in (4.3) is rewritten as an integral 



p u (k) = h 3 J dr' J drUl{r)U{r-r')p{r') 



-ikr 



(Bl) 
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where we used (4.1) and introduced an infinite mesh of Dirac delta functions 



m(r)= £ ^-rJ^E 



e -ik g m-r 



(B2) 



(We recall that fc 9 = 2n/h). Using the above representation of UI(r) and introducing in (Bl ) 
the Fourier series representation of the periodic charge density, 



(B3) 



fc'eK 



we recover the result (5.2) after straightforward simplifications. 



APPENDIX C: PROOF OF EQUIVALENCE BETWEEN EQS. (4.9) AND (4.6) 



Eq. (4.9) is equivalent to 



4S = ^E^( fc ) $ ( fc )' 



(Cl) 



fceK 



where $(fc) is the full Fourier transform (fc G K) of the back-interpolated potential 



mesh (4.8) 



$(fc) = / $(r)dr = h 3 dr e 



ikr 



V 



dr'm(r')U(r - r')<S> M (r'). 



(C2) 



We replace in this equation $M( r ') an d m(r') by their expressions (4.7) and (B2), and 
perform the integration over r': 



fc'eM m V 



(C3) 



The integration over r introduces a Kronecker symbol 5 k k / +kgrn . We get therefore the simple 
result 



(C4) 



where the function which is defined originally only for k e M, is now understood 

to be extended periodically to all K space. Notice that the inverse FFT does not introduce 



aliasing errors: the sum over m merely renders periodic. In accordance with (4.3) 
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and (4.5), we extend also Pm(^) an d G(k) periodically, with period 2ir/h. Using the above 



result and (4.5), the reciprocal energy rtClh can be expressed as 



4$ = ^£F(W*)m(*)G(fc) 

fceK 



P*(k + k g m)U(k + k g m)p M (k)G(k). 



2L 3 



(C5) 
(C6) 



fceM mez 3 



This may be compared with Eq. (4.6), i.e. 



P3M 2L 3 



^E^M( fc )PM(fc)G(fc). 



(C7) 



fceM 
fc^o 



Recalling (5.2) and the fact that U(k) is real, we see that both expressions are equivalent. 
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